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SUMMARY 


Current  nuclear  treaty  monitoring  interest  has  shifted  to  the  verification  of  a  comprehensive  test 
ban  treaty  [CTB7],  including  the  detection  of  the  first  test  of  potential  proliferalors.  Predictive  modeling 
capabilities  are  crucial,  not  just  desirable,  in  the  context  of  CTBT  and  proliferation  monitoring.  Synthetic 
seismograms  generated  with  the  linear-finite-difference  [LFD]  code  are  particularly  useful  for  regions 
where  earthquake  or  explosion  data  are  not  available.  However,  to  date  many  LFD  modeling  exercises 
are  still  limited  to  the  elastic  case  only. 

This  report  describes  a  simple  approach  to  include  anelastic  attenuation  linear-finite-difference  cal¬ 
culations.  The  algorithm  is  extremely  simple  to  implement,  and  it  gives  a  damping  effect  equivalent  to 
what  would  be  obtained  by  solving  the  viscoelastic  wave  equation.  Testings  of  this  algorithm  with  planar 
P,  Rg ,  and  Lg  waves  demonstrate  that  this  method  generates  a  frequency-independent  damping  over  a 
rather  broad  band,  which  is  equivalent  to  a  Q  increasing  linearly  with  frequency.  Any  Q  model  that  has 
a  form  other  than  Q  =  Qo  •  f^  can  be  simulated  easily  by  composing  LFD  results  generated  at  several 
separate  LFD  runs  correspoixiing  to  sampled  Q-f  pairs  of  the  desired  Q  structure.  Like  other  tech¬ 
niques,  this  algorithm  has  its  own  limitations.  One  shortcoming  is  that  the  performance  of  this  attenua¬ 
tion  operator  degrades  slightly  at  the  very  low  frequency.  Nevertheless,  the  shortcomings  of  this  algo¬ 
rithm  are  outweighed  by  the  simplicity.  Most  importantly,  this  procedure  preserves  causality. 

This  algorithm  has  been  applied  to  another  research  project  under  which  a  variety  of  mechanisms 
responsible  for  Lg  blockageAweakening  were  examined  and  compared.  With  the  help  of  this  new  algo¬ 
rithm,  it  is  rather  easy  to  equate  the  effects  due  to  small-scale  random  heterogeneities  (and  other  large- 
scale  structural  variations)  with  those  due  to  the  anelastic  ctttenuation.  The  results  are  summarized  in  a 
companion  report  "WsveguidB  offsets  of  Isrgo-scsio  structurol  vsriotion,  ono/ostic  sttonuotion,  ond  rsn- 
dom  heterogeneity  on  SV  Lg  propagation:  a  finite-difference  modeling  study”  (Report  PL-TR-96-2016). 


Received  for  publication  21  December  1995. 
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A  SIMPLE,  CAUSAL  METHOD  TO 
INCORPORATE  ANELASTIC  ATTENUATION  INTO 
FINITE-DIFFERENCE  CALCULATIONS 


1.  INTRODUCTION 

The  incorporation  of  intrinsic  attenuation  due  to  an  arbitrary  absorption  law  is  believed  to  be 
straightforward  in  frequency-domain  methods.  However,  so  far,  many  codes  used  in  waveform  synthesis 
(such  as  the  reflectivity  method  and  the  wavenumber-integration  method)  can  only  handle  a  frequency- 
independent  Q.  Bache  et  al.  (1981)  found  that  Lg  synthetics  made  for  models  with  frequency- 
independent  Q  do  not  attenuate  with  the  proper  frequency  dependence.  If  models  are  constructed  that 
match  Lg  amplitude-distance  relationship  at  1  Hz,  then  Lg  synthetics  would  attenuate  too  rapidly  at 
higher  frequencies. 

For  time-domain-based  numerical  methods,  it  has  been  very  difficult  to  add  the  intrinsic  attenuation 
because  the  anelastic  stress-strain  relation  has  the  form  of  a  convolution  integral,  which  is  intractable  in 
a  numerical  computation.  Vidale  and  Helmberger  (1988)  convolved  finite-difference  synthetics  with  a 
time-varying  operator  to  model  the  effect  of  anelastic  attenuation.  This  method  is  not  suitable  for  media 
with  spatially-varying  Q.  This  method  is  also  not  appropriate  when  there  are  significant  wave  conver¬ 
sions  (P  to  S  etc.)  in  media  where  Qp  is  not  equal  to  Qs-  The  first  successful  attempt  to  incorporate 
realistic  attenuation  laws  into  time-domain  methods  was  made  by  Day  and  Minster  (1984)  (see  also  Min¬ 
ster  et  al.,  1991)  based  on  the  method  of  Fade  approximation,  which  yields  an  expansion  of  the 
frequency-dependent  viscoelastic  modulus  into  a  rational  function.  Emmerich  and  Kom  (1987)  propose  a 
slightly  different  method  based  on  the  rheological  model  of  the  generalized  Maxwell  body,  which  has  a 
modulus  of  the  desired  rational  form.  The  major  inconvenience  of  these  approaches  is  that  they  demand 
a  continual  storage  of  five  or  more  time  steps  of  the  wave  field,  depending  on  the  accuracy  of  the 
approximation.  The  non-causal  nature  of  some  of  these  techniques  is  another  fundamental  dravirback,  at 
least  conceptually,  since  in  reality  the  anelastic  attenuation  of  the  Earth  should  act  in  a  causal  manner. 
That  is,  the  dissipation  of  energy  should  occur  as  soon  as  the  seismic  wave  arrives,  and  the  resulting 
displacement  at  the  current  time  step  should  not  be  dependent  on  that  of  the  future. 

With  all  these  considerations  in  mind,  a  different  procedure  is  developed  in  this  study  to  incor¬ 
porate  the  anelastic  attenuation.  It  turns  out  that  if  we  drop  the  ambitious  attempt  of  imposing  an  arbi¬ 
trary  (that  is,  user-defined)  frequency  dependence  on  the  Q  operator,  then  it  becomes  very  easy  to 
implement  a  causal,  phase-independent  damping  operator  which  is  quite  suitable  for  the  linear  finite- 
difference  (LFD)  calculation.  Several  researchers  have  demonstrated  that  the  performance  of  commonly 
used  absorbing  boundary  conditions  can  be  greatly  improved  if  a  viscous  damping  zone  is  added  to  the 
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grid  boundary  (Cetjan  et  al.,  1985;  Levander,  1985a).  The  damping  zone  simply  reduces  the  amplitude 
in  a  pointwise  manner.  There  is  no  reason  why  this  technique  cannot  be  exploited  to  model  the  anelastic 
attenuation. 

The  LFD  method  has  the  advantage  that  the  solution  contains  all  conversions  and  all  orders  of 
multiple  scattering.  It  permits  examinations  of  fairly  general  models  with  arbitrary  complex  variations  in 
material  properties  and  free-surface  geometry.  Furthermore,  it  does  not  require  many  assumptions  com¬ 
monly  invoked  in  other  theoretical  approaches.  The  basic  limitations  to  the  LFD  method  are  the  compu¬ 
tational  cost  and  memory  requirements.  These  constrain  the  size  of  the  grid  and  the  number  of  time 
steps  that  can  be  calculated  in  a  reasonable  time.  Despite  the  limitations,  LFD  is  still  one  of  the  most 
powerful  and  popular  modeling  tools  in  generating  realistic  synthetic  seismograms,  and  it  is  extremely 
desirable  to  have  a  convenient  means  to  incorporate  anelastic  attenuation  into  LFD  calculations.  In  this 
report,  a  causal  procedure  for  this  purpose  is  derived  and  discussed,  followed  by  three  sets  of  LFD 
experiments.  Different  seismic  waves  are  used  as  the  initial  pulse  in  these  examples  to  demonstrate  the 
effectiveness  of  this  algorithm. 
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2.  METHOD 


Consider  the  simplest  isotropic  homogeneous  medium  in  which  the  amplitude  of  seismic  waves 
decays  exponentially  with  travel  distance: 

A(f,A)sAo(f)  G*Ao(f)  e‘"  ""  ,  [1] 


yi^ere  y  »  JLlL^  U  js  the  group  velocity,  and  A  is  the  distance  traveled.  In  LFD  calculations,  A  is  taken 
to  be  the  distance  that  the  seismic  wave  would  travel  within  one  temporal  step  of  the  LFD  iteration,  that 

n  •  1  ♦  dt 

is,  A  =  U  •  dt.  Thus  G  can  also  be  written  as  e  °  -  There  exists  a  constant  such  that 

Gae“^  "^  =  1 -tr  y  A.  [2] 

It  Implies  that  G  =  1  -  ^  -  and  hence  Q  =  Qo  -  f  if  we  define  Qo  to  be  If  the  damp¬ 

ing  term  G  Is  a  function  of  the  grid  coordinate  only  and  invariant  with  frequency,  then  we  would  have  a 
Q  increasng  linearly  with  frequency.  If,  however,  a  frequency-dependent  G  is  used  in  separate  LFD 
simulations,  then  combining  the  band-limited  LFD  results  would  produce  the  solution  for  that  particular 
frequency-dependent  Q  model.  Here  the  parameter  q  is  a  function  of  y  or  G,  and  the  means  of  deter¬ 
mining  Ti  will  be  discussed  later.  In  practice,  however,  the  users  only  need  to  specify  a  multiplicative 
constant  G  slightly  less  than  1  for  each  grid  point  These  localized  damping  factors  are  used  to  modify 
the  displacement  field  pointwise  at  each  iteration  step.  The  decay  rate  (y)  and  the  quairty  factor  (Q)  can 
be  determined  later  after  the  finite-difference  calculation  is  done.  A  possible  drawback  of  this  approach 
is  that  if  a  specific  frequency-dependent  Q  model  is  desired,  then  several  separate  LFD  simulations 
need  to  be  carried  out  for  each  frequency-Q  pair  of  interest,  as  discussed  above.  Nevertheless,  this 
possible  shortcoming  of  this  approach  is  outweighed  by  its  simplicity.  More  importantly,  this  procedure 
preserves  the  causality.  Another  characteristic  of  this  approach  is  that,  given  a  damping  factor  G  in  a  P- 

SV  LFD  calculation,  Qo  =  would  be  applicable  to  both  P  and  S  phases.  Thus  Qp  and  Qs 

should  be  about  the  same. 


So  far  we  have  derived  several  necessary  conditions  for  an  anelastic  attenuation  model,  based  on 
the  desired  exponential  decay  of  seismic  amplitude.  In  the  following,  we  shall  take  a  schematic  view  of 
this  proposed  algorithm.  Consider  the  heterogeneous  acoustic  wave  equation  in  the  nondissipative 
medium: 


a^p 

at^ 


c2{x,  z)[  ^^2 


a^p, 


where  P  is  the  acoustic  displacement  potential,  and  c(x,  z)  is  the  acoustic  velocity  at  the  node  (x. 


[3] 

z)- 
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LFD  iteratively  solves  for  the  current  pressure  at  (x,  z)  using  the  past  pressure  within  a  spatially-limited 
region  surrounding  the  grid  (x,  z).  The  actual  size  of  the  temporal  and  spatial  buffers  required  depends 
on  the  order  of  the  LFD  scheme.  At  each  time  step,  the  resulting  unattenuated  pressure  P  is  then  multi¬ 
plied  by  the  damping  G  to  get  the  dissipated  f*.  Once  the  whole  pressure  field  is  updated  and  dissi¬ 
pated  (grid  by  grid)  with  the  corresponding  spatial  damping  factors,  the  standard  LFD  iteration  restarts 
for  the  next  pressure  field  without  damping,  and  then  damping  factors  are  applied  again.  The  procedure 
for  the  P-SV  LFD  calculation  is  exactly  the  same.  It  should  be  obvious  that  the  extra  calculation  and 
memory  required  in  this  procedure  are  minimum.  We  can  envision  the  pointwise  damping  factor  as  a 
degenerate  digital  filter  which  has  only  one  point  of  temporal  span  (that  is,  memory).  It  can  be  regarded 
as  the  limiting  case  of  the  convolutional  integral  of  many  attenuation  operators  typically  used  in  the 
frequency-domain  approach. 

The  attenuated  pressure,  =  P  ■  G,  can  also  be  written  as 

P  =  P  -  ^  ^  •  P  .  [4] 

That  is,  P  can  be  obtained  by  adjusting  the  undissipated  P  a  little.  The  coefficient  of  pressure  loss, 
" ■  is  very  similar  to  the  coefficient  of  the  friction  term  which  Levander  (1985a)  (and  Frankel  and 

Wennerberg,  1987)  used  in  the  telegraphy  equation.  The  difference  is  that  we  have  added  the  t)  term 
here  to  obtain  the  correct  relationship  between  Q  and  the  damping  effect. 
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3.  EXAMPLE  1:  P-WAVE  EXPERIMENT 


The  first  example  includes  four  separate  LFD  runs  to  test  the  propagation  of  a  normally  incident  P 
wave  through  one  elastic  and  three  anelastic  models.  The  whole  space  has  compressional  and  shear 
velocities  of  6.20  and  3.58  km/sec,  respectively.  The  LFD  grid  has  500  nodes  in  the  Z  direction,  and  an 
absorbing  boundary  condition  is  applied  to  the  top  of  the  grid.  The  spatial  and  temporal  spadngs  are 
250  meters  and  25  ms,  respectively.  In  each  test,  the  initial  incident  wave  is  a  planar  pulse  with  a  shape 
of  the  derivative  of  Gaussian  curve.  The  pulse  is  32  grids  long  and  centered  at  1.5  Hz.  For  the  three 
anelastic  models,  the  damping  factors  are  arbitrarily  set  to  (A)  0.995,  (B)  0.999,  and  (C)  0.9995,  respec¬ 
tively.  Figure  1  shows  the  vertical-component  record  sections  of  both  the  elastic  model  (top)  and  the 
anelastic  model  A  with  a  damping  factor  of  0.995  (bottom),  plotted  in  the  same  scale.  In  this  example 
(and  other  examples  given  in  latter  sections),  the  geometrical  spreading  is  not  present  in  the  LFD  solu¬ 
tion,  and  hence  the  standard  one-station  spectral  analysis  can  be  used  to  measure  the  amplitude  decay 
and  Q  without  the  need  to  correct  for  the  geometrical  spreading  first.  A  signal  window  of  25.6  seconds 
(1024  points)  is  used  to  compute  the  spectral  ratio  and  the  quality  factor  (Figures  2  and  3  show  exam¬ 
ples).  The  elastic  model  gives  a  y  less  than  0.0001 ,  which  can  be  regarded  as  a  measure  of  the  numer¬ 
ical  accuracy  of  the  2nd-order  explicit  LFD  scheme  used  here.  Higher-order  LFD  schemes  should  yield  a 
Y  even  smaller  for  the  elastic  case.  Rgure  4  shows  the  spectral  ratio  of  trace  No.  1  to  trace  No.  16  of 
the  model  A.  These  two  sensors  are  75  km  apart.  It  is  interesting  to  note  that  actually  the  spectral  ratio 
(and  hence  the  resulting  y)  is  nearly  a  constant  over  a  rather  broad  band.  The  y  of  this  model,  meas¬ 
ured  at  a  suite  of  sensors  evenly  spread  over  a  wide  range,  is  consistently  0.018.  The  y  values  of 
models  B  and  C  are  0.003  and  0.002,  respectively.  Table  1  below  lists  the  y  and  quality  factors  meas¬ 
ured  at  fourteen  sensors.  In  each  case,  the  resulting  Q  is  clearly  a  linear  function  of  frequency. 

Futterman  (1962)  showed  that  if  a  linear  theory  of  wave  propagation  is  assumed,  then  the  pres¬ 
ence  of  absorption  is  a  necessary  and  sufficient  condition  for  the  presence  of  dispersion.  He  gave  an 
excellent  discussion  of  absorption-dispersion  pairs.  Using  the  condition  of  causality,  he  derived  a 
Kramers-Krong  relation  that  allows  one  to  calculate  the  dispersion  for  a  given  absorption  function.  (See 
also  Ganley  1981;  Liu  et  al.,  1976;  Jacobson,  1987.)  For  frequency-independent  Q  models,  this  often 
leads  to  a  noticeable  broadening  in  the  signal  pulse  as  the  seismic  pulse  travels  to  longer  distances.  In 
the  special  case  of  frequenry-independent  damping  (y),  or  equivalently,  Q  =  Qo  •  f  (such  as  the  damp¬ 
ing  we  used  here),  there  should  be  no  significant  broadening.  The  dispersion  is  present,  as  a  combined 
result  of  the  causality  and  the  grid  dispersion.  This  is  illustrated  in  Figure  1 . 
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Model 

ID.G 


Table  1.  P-wave  Experiment  with  Three  Anelastic  Models 


2.5  Hz 


Model 

Sensor 

1  Hz 

1.5  Hz 

2  Hz 

2.5  Hz 

ID.G 

Range  (km) 

Qo 

Y 

Q 

Y 

Q 

Y 

Q 

Y 

A.  0.995 

75 

28 

42 

0.018 

59 

0.017 

75 

0.017 

A,  0.995 

70 

28 

■TibB 

42 

0.018 

59 

0.017 

75 

0.017 

A,  0.995 

65 

28 

0.018 

42 

0.018 

59 

75 

0.017 

A,  0.995 

60 

28 

0.018 

42 

0.018 

59 

75 

0.017 

A.  0.995 

55 

28 

0.018 

42 

0.018 

59 

75 

0.017 

A,  0.995 

50 

28 

0.018 

42 

0.018 

59 

0.017 

75 

0.017 

A,  0.995 

45 

28 

0.018 

42 

0.018 

59 

0.017 

75 

0.017 

A,  0.995 

40 

28 

0.018 

42 

0.018 

59 

0.017 

75 

0.017 

A,  0.995 

35 

28 

0.018 

42 

0.018 

59 

0.017 

75 

0.017 

A.  0.995 

30 

28 

0.018 

42 

0.018 

59 

0.017 

75 

0.017 

A.  0.995 

25 

28 

0.018 

42 

0.018 

59 

0.017 

75 

0.017 

A,  0.995 

20 

28 

0.018 

42 

0.018 

59 

0.017 

75 

0.017 

A,  0.995 

15 

28 

0.018 

42 

0.018 

59 

75 

0.017 

A,  0.995 

10 

28 

0.018 

42 

0.018 

59 

0.017 

75 

0.017 

B,  0.999 

75 

156 

0.003 

228 

0.003 

303 

0.003 

369 

B,  0.999 

70 

155 

0.003 

228 

0.003 

304 

0.003 

388 

B,  0.999 

65 

154 

0.003 

229 

0.003 

305 

0.003 

387 

0.003 

B.  0.999 

60 

154 

0.003 

229 

0.003 

306 

0.003 

378 

0.003 

B,  0.999 

55 

154 

0.003 

229 

0.003 

307 

0.003 

380 

0.003 

B,  0.999 

50 

154 

0.003 

229 

0.003 

308 

0.003 

386 

0.003 

B.  0.999 

45 

154 

0.003 

229 

0.003 

308 

0.003 

384 

0.003 

B,  0.999 

40 

154 

0.003 

229 

0.003 

308 

0.003 

379 

0.003 

B,  0.999 

35 

154 

0.003 

229 

0.003 

308 

0.003 

381 

0.003 

B,  0.999 

30 

154 

0.003 

229 

0.003 

307 

0.003 

390 

0.003 

B,  0.999 

25 

154 

0.003 

229 

0.003 

306 

0.003 

383 

0.003 

B,  0.999 

20 

154 

0.003 

228 

0.003 

304 

0.003 

375 

0.003 

B,  0.999 

0.003 

225 

0.003 

303 

0.003 

383 

0.003 

B,  0.999 

0.003 

228 

0.003 

305 

0.003 

380 

0.003 

C.  0.9995 

75 

323 

0.002 

466 

0.002 

608 

0.002 

708 

C,  0.9995 

70 

318 

0.002 

464 

0.002 

0.002 

787 

C.  0.9995 

65 

316 

0.002 

464 

0.002 

610 

0.002 

792 

C,  0.9995 

60 

317 

0.002 

465 

0.002 

613 

0.002 

755 

C.  0.9995 

55 

317 

0.002 

463 

0.002 

616 

0.002 

753 

0.002 

C.  0.9995 

50 

318 

0.002 

464 

0.002 

620 

0.002 

778 

0.002 

C,  0.9995 

45 

316 

0.002 

464 

0.002 

625 

0.002 

786 

0.002 

C.  0.9995 

40 

315 

0.002 

460 

0.002 

626 

0.002 

760 

0.002 

C,  0.9995 

35 

317 

0.002 

462 

0.002 

631 

0.002 

751 

0.002 

C,  0.9995 

30 

314 

0.002 

460 

0.002 

629 

0.002 

799 

0.002 

C,  0.9995 

25 

317 

0.002 

467 

0.002 

629 

0.002 

801 

0.002 

C.  0.9995 

20 

313 

0.002 

459 

0.002 

615 

0.002 

745 

0.002 

C,  0.9995 

15 

303 

0.002 

442 

0.002 

601 

0.002 

772 

0.002 

C,  0.9995 

10 

313 

0.002 

457 

0.002 

613 

0.002 

764 

0.002 

It  was  shown  in  the  previous  section  that  Qo  =  -^G^’  ®  damping  factor  G  in  a  P-SV 

LFD  calculation.  Since  this  Qo  formula  is  applicable  to  both  P  and  S  phases,  Qp  and  Os  should  be 
about  the  same.  This  can  easily  be  tested  by  using  identical  grid  parameters  in  the  LFD  experiments  but 
with  different  initial  pulse  types,  as  shown  in  Table  2  below. 


Table  2.  S-wave  Experiment  with  Three  Anelastic  Models 


Model 

Sensor 

1  Hz 

1.5  Hz 

2  Hz 

2.5  Hz 

ID,  G 

Range  (km) 

Y 

o 

Y 

Q 

Y 

Q 

Y 

A.  0.995 

200 

27 

0.033 

El 

0.033 

56 

0.031 

75 

0.029 

A,  0.995 

160 

27 

0.033 

n 

0.033 

56 

0.031 

71 

0.031 

A,  0.995 

120 

26 

0.034 

40 

0.033 

57 

0.031 

70 

0.031 

A,  0.995 

80 

26 

0.034 

40 

0.033 

57 

0.031 

71 

0.031 

A,  0.995 

40 

26 

0.035 

40 

0.033 

57 

0.031 

71 

0.031 

B.  0.999 

200 

146 

0.006 

0.006 

299 

0.006 

362 

0.006 

B,  0.999 

160 

146 

0.006 

0.006 

297 

0.006 

363 

0.006 

B.  0.999 

120 

146 

0.006 

223 

0.006 

297 

0.006 

363 

0.006 

B,  0.999 

80 

146 

0.006 

223 

0.006 

297 

0.006 

363 

0.006 

B,  0.999 

40 

145 

0.006 

0.006 

297 

0.006 

362 

0.006 

C,  0.9995 

200 

303 

0.003 

453 

0.003 

611 

0.003 

734 

0.003 

C,  0.9995 

160 

304 

0.003 

451 

0.003 

598 

0.003 

726 

0.003 

C,  0.9995 

120 

304 

0.003 

451 

0.003 

598 

0.003 

726 

0.003 

C.  0.9995 

80 

305 

0.003 

451 

0.003 

598 

0.003 

727 

0.003 

C,  0.9995 

40 

304 

0.003 

451 

0.003 

598 

0.003 

727 

0.003 

Based  on  the  results  in  Tables  1  and  2.  we  can  solve  for  n  in  each  case  (Table  3).  It  should  be 
emphasized  (again)  however,  that  the  parameter  ii  is  introduced  merely  for  convenience  in  deriving  the 
algorithm.  In  practice,  it  is  not  necessary  to  determine  the  value  of  ti  in  advance  to  incorporate  the 
damping  effect  into  the  LFD  calculation. 


Table  3.  Damping  Factor  and  Associated  Qq  and  t| 


Damping 

P-wave  Result 

S-wave  Result 

G 

Qo(P) 

tl 

Qo(S) 

n 

0.9950 

28 

5.6/71 

27 

0.9990 

154 

6.1/71 

146 

0.9995 

317 

6.3/n 

303 

6.1/n 
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1 

1 

1 

1 

1 

1 

1 

1 

1 

1 

1 

1 

1 

1 

1 

1 

1 

1 

1 

1 


1 

1 

1 

1 

1 

1 

1 

1 

1 

1 

1 

1 

1 

1 

1 

1 

1 

1 

1 

1 


-5km 

0.439E-^O0 

-10km 

0.493E+00 

-15km 

0.492E4O0 

-20km 

0.494E-FO0 

-25km 

0.493E+O0 

-30km 

0.493E+00 

-35km 

0.493E+00 

-40km 

0.492E400 

•45km 

0.492E+00 

-50km 

0.488E-K)0 

-55km 

0.488E400 

-60km 

0.483E-F00 

-65km 

0.482E400 

-70km 

0.476E+00 

-75km 

0.473E*»00 

-80km 

0.465E4O0 

-85km 

0.461  E-»«0 

-90km 

0.451  E+00 

-95km 

0.447E400 

-100km 

0.436E-fO0 

Second 


P  wave,  model  no_Q,  dx=dz=250m,  dt=:25ms 

VERT,  same  scale,  decimated  by  3 


-5km 

0.571  E-01 

-10km 

0.672E-01 

-15km 

0.756E-01 

-20km 

0.834E-01 

-25km 

0.870E-01 

-30km 

0.972E-01 

-35km 

0.113E-fO0 

-40km 

0.125E+00 

-45km 

0.133E400 

-50km 

0-143E+00 

-55km 

0.170E-^00 

-60km 

0.192E+00 

-65km 

0.205E400 

-70km 

0.216E-FO0 

-75km 

0.262E400 

-80  km 

0.300E400 

-85km 

0.318E400 

-90km 

0.343E400 

-95km 

0.433E+00 

-100km 

0.475E-fO0 

Second 


P  wave,  model  Q.995,  dx=dz=250m,  dt=25ms 

VERT,  same  scale,  dectmaied  by  3 


Figure  1.  Vertical-component  record  sections  of  normally-incident  planar  P  waves  in  the  elastic  (top)  and  aneiastic  (bot¬ 
tom)  whole-space  models,  respectively.  Sensors  are  5  km  apart 
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FREQUENCY  (Hz) 


Specalc+Waveshaper 

version  9502,  JB+TWM+RJ 

Window  1 024 
Signal  1  1024 

15.00%  Cosine  Taper 
Taper  Both  Ends 

41  R  Smoothing 
No  Decimation 

Deconvolution  Flag=0. 

Waveshaper  Flag=0. 

Reject:  <  O.OOHz 
Reject:  >1  O.OOHz 
PassBand=  0.01.10.00Hz 

Final  Taper=1 0.00%, 10.00% 

Model  P.Q.995 
VERT.16.  Pwave 


1 2.79  sec  0.3004  Input  Trace  1 2.79  sec  0.0003  Output  T race 


Figure  2.  Amplitude  spectrum  of  the  vertical-component  LFD  synthetic  displacement  seismogram  recorded  at  sensor 
No.  16  of  the  anelastic  model  (A)  which  has  a  damping  factor  of  0.995. 
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Specalc+Waveshaper 

version  9502,  JB+TWM+RJ 

Window  1 024 
Signal  481  1024 

15.00%  Cosine  Taper 
Taper  Both  Ends 

41  R  Smoothing 
No  Decimation 

Deconvolution  Flag=0. 

Waveshaper  Flag=0. 

Reject:  <  O.OOHz 
Reject:  >1  O.OOHz 
Pass  Band=  0.01,10.00Hz 

Final  Taper=  10.00%,  10.00% 

Model  P.Q.995 
VERT.1 ,  P  wave 


1 .00  2.00  3.00 

FREQUENCY  (Hz) 


12.79  sec  0.0571  Input  Trace 


12.79  sec  0.0001  Output  Trace 


Figure  3.  Same  as  Figure  2  except  for  sensor  No.  1.  which  is  75  km  beyond  sensor  No.  16. 
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TGAL  RATIO  vers  1.5  (JB+TWM+RAW+RJ) 

Thu  Nov  9  08:07:24  1995 

Model  P.Q.995.  ratio:  No.1  /  No.1 6, 75  km  apart 

NOISE  POWER  NOT  SUBTRACTED  FROM  SIGNAL  POWER 

S/N  POWER  THRESHOLD  =  2.0 

NO  SMOOTHING 


O 


o 


FREQ  BAND  SLOPE  +/-  Y INTERC  AV INTEG  AV  RATIO  SD 
0.75-1.25  2.175e-02  1.993e-03  2.441  e-01  -5.906e-01  -5.903e-01  3.498e-03 

1.75-2.25  1.622e-02  1.483e-03  2.563e-01  -5.588e-01  -5.589e-01  2.599e-03 

1.00-2.00  3.685e-02  1.608e-03  2.341e-01  -5.754e-01  -5.752e-01  1.128e-02 

2.00-3.00  7.322e-03  2.015e-03  2.682e-01  -5.532e-01  -5.532e-01  3.577e-03 

0.80-2.40  3.114e-02  8.365e-04  2.387e-01  -5.723e-01  -5.723e-01  1.479e-02 

Figure  4.  Spectral  amplitude  ratio  (shown  in  log  scale)  of  synthetic  seismograms  No.  1  to  No.  16  (cf. 
Figures  2  and  3)  of  the  anelastic  model  (A)  which  has  a  damping  factor  of  0.995.  The  amplitude  ratio  is 
nearly  a  constant  over  a  rather  broad  band,  which  corresponds  to  a  Q  increasing  linearly  with  frequency 
(see  Table  1). 
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4.  EXAMPLE  2:  Lg-WAVE  EXPERIMENT 


In  the  second  example  we  test  the  applicability  of  the  proposed  method  to  Lg  waves  with  a 
single-layer  crustal  model.  The  pure  Lg  wave  packet  is  injected  into  a  stratified  portion  of  the  LFD  grid 
as  the  reference  initial  condition  to  trigger  all  LFD  calculations.  An  obvious  advantage  of  this  approach  is 
that  the  effect  due  to  different  types  of  heterogeneity  on  Lg  phase  can  be  easily  isolated  and  evaluated. 
The  algorithm  generating  the  incident  Lg  wave  packet  can  be  found  in  Jih  (1995)  and  it  is  analogous  to 
the  one  Boore  (1970)  developed  for  the  fundamental-mode  Love  wave  packet.  The  same  approach  has 
also  been  extensively  used  in  elastic  (P-SV)  LFD  simulations  of  Rg  propagation/scattering  problems. 
Among  many  others  using  the  approach  are:  Levander  (1985b).  Toksoz  et  al.  (1986),  McLaughlin  and 
Jih  (1986, 1987),  and  Jih  (1993b,  1996). 

On  real  seismograms,  the  Lg  phase  often  lacks  a  dear  onset,  but  it  does  have  a  well-ddined 
amplitude  maximum  with  a  group  velodty  around  3.5  km/sec.  The  Lg  waves  are  basically  the  interfer¬ 
ence  of  multiply  reflected  S  waves  boundng  back  and  forth  between  the  free  surface  and  the  Moho. 
We  can  envision  that  a  system  of  planar  S  waves  is  set  off  at  the  Moho  at  equal  delay  with  the  same 
post-critical  indinadon,  while  another  system  of  planar  S  waves  is  set  off  at  the  free  surface  in  a  sym¬ 
metric  manner,  which  corresponds  to  the  bundle  of  reflected  waves.  If  these  two  systems  of  waves  (or 
rays)  are  properly  confined  in  the  same  stratified  region  of  the  crust,  they  will  establish  repetitive 
reflections  (at  both  the  free  surface  and  the  Moho)  with  constructive  interferences.  By  adjusting  the  time 
delay  between  the  consecutive  5  wavefronts,  we  can  obtain  a  denser  (or  coarser)  interference  pattern. 
The  Lg  waves  can  also  be  described  as  the  superposition  of  many  higher-mode  surface  waves  that 
interfere  to  give  the  complex  observed  waveforms.  These  modes  correspond  to  waves  trapped  in  the 
crustal  wave  guide. 

The  pure  Lg  wave  packet  generated  with  the  procedure  described  in  Jih  (1995, 1996)  exhibits  all 
these  expeded  features.  Figure  5  gives  the  vertical-component  snapshots  of  Lg  wave  propagation  in  an 
anelastic  single-layer  crustal  model  taken  at  a  temporal  spadng  of  10  seconds.  The  homogeneous  crust 
is  30  km  thick  with  P-  and  S-wave  velocities  of  6.2  and  3.58  km/sec,  respectively.  The  damping  factor 
in  the  crust  is  set  to  0.999.  Since  no  scattering  mechanism  is  present,  the  checkerboard-like  pattern  due 
to  the  constructive  interference  of  repeatedly  reflected  5  waves  trapped  in  the  crust  is  retained  at  all 
times  (Rgure  5).  This  chessboard-like  pattern  undoubtedly  indicates  that  the  two  commonly  quoted 
interpretations  of  Lg  waves,  either  as  multiply  reflected  S  waves  or  as  higher-mode  surface  waves,  are 
indeed  adequate.  Note  that  the  computational  grid  was  shifted  at  15  seconds  (not  shown)  to  stretch  the 
lateral  span  of  the  grid.  The  “marching  grid”  technique  is  described  in  Jih  (1993a). 

Figure  6  displays  the  record  section  of  the  vertical-component  synthetics  of  this  anelastic  model. 
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For  comparison,  the  record  section  corresponding  to  the  elastic  case  is  also  included.  We  fmd  that 
several  parameters,  Including  the  group  velocity  and  the  attenuation,  can  be  more  reliably  measured 
with  synthetics  recorded  at  greater  depths.  Figure  7  shows  the  record  sections  at  a  depth  of  10  km. 
Based  on  these  seismic  sections,  the  group  velocity  of  Lg  wave  packet  is  determined  to  be  3.33 
km/sec,  which  is  appropriate  for  areas  like  the  western  U.S. 

Table  4  lists  the  quality  factors  and  amplitude  ratios  measured  at  a  suite  of  sensors  spread  from 
130  to  180  km  at  two  depths.  For  comparison,  the  results  for  two  less  attenuative  models  with  damping 
factors  0.9995  and  0.9998,  respectively,  are  also  included.  For  each  model,  the  quality  factors  are 
essentially  constant  across  the  array  at  the  same  depth.  However,  Q  values  measured  at  the  free  sur¬ 
face  are  systematically  smaller  than  those  measured  at  10  km.  Also,  the  y  values  at  2  Hz  are  slightly 
smaller  than  those  of  1  Hz.  Nevertheless,  it  is  dear  that  the  Q  we  modeled  can  still  be  represented 
approximately  as  a  linear  function  of  frequency.  Figures  8  through  10  show  a  set  of  amplitude  spectra 
and  the  corresponding  amplitude  ratio  for  the  anelastic  model  with  a  damping  factor  of  0.999.  Over  the 
frequency  band  0.8-2.4  Hz,  the  ratio  is  rather  flat 
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Table  4.  Lg  Propagation  in  Anelastic  Crustal  Models 


Damper 

Sensor  Location  (km) 

Results  at  1  Hz 

Results  at  2  Hz 

G 

Range 

Depth 

A/Ao 

Qo 

Yd  Hz) 

A/Ao 

0(2  Hz) 

y(2Hz) 

0.9990 

130 

0 

0.286 

100 

0.0096 

239 

0.9990 

140 

0 

0.262 

100 

0.0096 

243 

0.9990 

150 

0 

0.242 

101 

0.0095 

232 

0.0083 

0.9990 

160 

0 

0.230 

104 

0.0092 

0.301 

259 

0.0074 

0.9990 

170 

0 

0.227 

110 

0.0087 

0.300 

276 

0.0070 

0.9990 

180 

0 

0.201 

108 

0.0089 

0.268 

267 

0.0072 

0.9990 

130 

10 

0.352 

120 

0.0080 

0.377 

251 

0.9990 

140 

10 

0.316 

116 

0.0082 

0.359 

261 

0.9990 

150 

10 

0.298 

119 

0.0081 

0.333 

261 

0.0073 

0.9990 

160 

10 

0.284 

122 

0.0079 

0.309 

262 

0.0073 

0.9990 

170 

10 

0.274 

126 

0.0076 

0.289 

262 

0.0073 

0.9990 

180 

10 

0.252 

125 

0.0077 

0.259 

253 

0.0076 

0.9995 

130 

0 

0.489 

174 

0.0055 

0.555 

420 

0.0046 

0.9995 

140 

0 

0.460 

0.0055 

0.538 

433 

0.0044 

0.9995 

150 

0 

0.441 

^^9 

0.0055 

0.499 

403 

0.0048 

0.9995 

160 

0 

0.440 

187 

0.0051 

0.526 

486 

0.0039 

0.9995 

170 

0 

0.438 

198 

0.0048 

0.534 

537 

0.0036 

0.9995 

180 

0 

0.410 

194 

0.0049 

0.505 

516 

0.0037 

0.9995 

130 

10 

0.589 

235 

0.0041 

0.601 

474 

0.0040 

0.9995 

140 

10 

0.551 

225 

0.0043 

0.591 

502 

0.0038 

0.9995 

150 

10 

0.535 

230 

0.0042 

0.566 

500 

0.0038 

0.9995 

160 

10 

0.526 

239 

0.0040 

0.545 

507 

0.0038 

0.9995 

170 

10 

0.525 

253 

0.0038 

0.523 

496 

0.0039 

0.9995 

180 

10 

0.495 

246 

0.0039 

0.482 

462 

0.0041 

0.9998 

130 

0 

0.669 

310 

0.0031 

0.730 

772 

0.0025 

0.9998 

140 

0 

0.644 

305 

0.0031 

0.719 

808 

0.0024 

0.9998 

150 

0 

0.625 

306 

0.0031 

0.682 

714 

0.0027 

0.9998 

160 

0 

0.642 

347 

0.0028 

0.736 

1034 

0.0019 

0.9998 

170 

0 

0.649 

378 

0.0025 

0.759 

1272 

0.0015 

0.9998 

180 

0 

0.620 

362 

0.0026 

0.729 

1135 

0.0017 

0.9998 

130 

10 

0.799 

556 

0.0017 

0.793 

992 

0.0019 

0.9998 

140 

10 

0.764 

500 

0.0019 

0.794 

1106 

0.0017 

0.9998 

150 

10 

0.755 

512 

0.0019 

0.776 

1095 

0.0018 

0.9998 

160 

10 

0.757 

552 

0.0017 

0.762 

1135 

0.0017 

0.9998 

170 

10 

0.768 

617 

0.0016 

0.741 

1053 

0.0018 

0.9998 

180 

10 

0.741 

574 

0.0017 

0.698 

910 

0.0021 
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FDstep  0 

VERT  OiXXtsec,  1  800  1  500 
XTIcMark=  10  km;  Z  Tic  Mark  =  10km 
-0^6E+00,  0^6E+00 


FDstep  400 

VERTIOilOOsec,  1  800  1  500 
XTlcMark=  10  km;  Z  Tic  Mark  =  10km 
-0J2BE*00,  0^*00 


FDstep  800:  grid  shifted 
VERT  20i)00sec,  1  800  1  500 

XTIcMarks  10  km;  Z  Tic  Mark  =  10km 
-024E-ID0,  O^E^OO 


FDstep  1200 

VERT30.000SOC,  1  800  1  500 
XTIcMarkr  10km;  ZTlcMarks  10km 
-O.ISE-hOO,  0.19E+00 


FDstep  1600 

VERT40.0008ec,  1  800  1  500 

XTicMaiks  10  km;  Z  Tic  Mark  =  10km 
-0.14E-M)0,  0.14E-ID0 


FDstep  2000 

VERT50J)00sec,  1  800  1  500 
XTlcMark=  10  km;  Z  Tic  Marks  10km 
>0.11E4>00,  0.11E400 


FDstep  2400 

VERT  eOUXXtooc,  1  800  1  500 
XTicMarks  10  km;  Z  Tic  Marks  10km 
*0^7E-01»  O^E-01 


LFD  Simulation  of  Lg  Propagation:  Model  Q.999 

Figure  5.  The  vertical-component  snapshots  of  Lg  visNe  propagation  in  an  1 -layer  anelastic  crustal 
model.  The  maximum  amplitude  of  the  wavefield  (shown  at  the  right)  drops  from  0.360  to  0.094  after  trav¬ 
eling  a  distance  of  200  km.  The  chesstx)ard-like  pattern  of  the  Lg  wave  packet,  which  is  due  to  the  con¬ 
structive  interference  of  multiply  reflected  planar  S  waves  in  the  crust,  is  retained,  however. 
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Rgure  6.  Vertical-component  record  sections  of  planar  Lg  wave  propagation  recorded  at  surface  sensors  (10  km  apart) 
for  the  1 -layer  elastic  (top)  and  anelastic  (bottom)  crustal  models,  respectively. 
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Figure  7.  Vertical-component  record  sections  of  planar  Lg  waves  recorded  at  sensors  10  km  deep  for  the  1 -layer  elas¬ 
tic  (top)  and  anelastic  (bottom)  crustal  models,  respectively.  Parameters  measured  off  synthetics  at  this  depth,  such  as 
the  peak  amplitude  and  the  group  velocity,  are  more  reliable  than  those  measured  at  the  free  surface. 
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version  9502,  JB+TWM+RJ 

Window  1 024 
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1 0.00%  Cosine  Taper 
Taper  Both  Ends 

41  R  Smoothing 
Decimated  By  2 

Deconvolution  Flag=0. 

Waveshaper  Flag=0. 

Reject:  <  O.OOHz 
Reject:  >1  O.OOHz 
Pass  Band=  0.01,10.00Hz 

Final  Taper=1 0.00%,  1 0.00% 

Model  Q.999 
VERT.34,  Lg 


12.79  sec  0.2626  Input  Trace  12.79  sec  0.0003  Output  Trace 


Rgure  8.  Amplitude  spectrum  of  the  vertical-component  displacement  recorded  at  sensor  No.  34  of  the  anelastic  model. 
The  sensor  is  located  near  the  left  portion  of  the  grid,  at  a  depth  of  10  km. 
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Figure  9.  Amplitude  spectrum  of  the  vertical-component  displacement  recorded  at  sensor  No.  52  of  the  anelastic  model. 
The  sensor  is  located  near  the  right  portion  of  the  grid,  at  a  depth  of  1 0  km. 
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TGAL  RATIO  vers  1.5  (JB+TWM+RAW+RJ) 

Thu  Nov  9  08:20:03  1995 
Model  Q.999,  ratio:  No.52  /  No.34 

NOISE  POWER  NOT  SUBTRACTED  FROM  SIGNAL  POWER 
S/N  POWER  THRESHOLD  =  2.0 
7  POINT  SMOOTHING 


O 


FREQ  BAND  SLOPE  +/-  Y INTERC  AV INTEG  AV  RATIO  SD 
0.75-1.25  8.689e-04  3.475e-04  2.515e-01  -5.985e-01  -5.985e-01  1.993e-04 

1.75-2.25  6.625e-03  1.528e-02  2.509e-01  -5.872e-01  -5.872e-01  5.559e-03 

1.00-2.00  1.479e-02  2.644e-03  2.427e-01  -5.927e-01  -5.930e-01  5.356e-03 

2.00-3.00  5.150e-02  1.182e-02  2.051e-01  -5.593e-01  -5.593e-01  2.034e-02 

0.80-2.40  1.554e-02  2.074e-03  2.424e-01  -5.906e-01  -5.906e-01  8.369e-03 


Rgure  10.  Lg  amplitude  ratio  of  synthetic  seismograms  No.  52  to  No.  34  of  the  anelastic  model.  The 
sensors  are  180  km  apart  (see  Figures  7,  8.  and  9).  The  amplitude  ratio  is  nearly  a  constant  of  0.24  over 
a  rather  broad  band. 
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5.  EXAMPLE  3:  Rg-WAVE  EXPERIMENT 


Recent  observational  studies  indicate  that  Rg  could  play  an  important  role  in  Lg  excitation  for 
Yucca  Flat  explosions  (Patton  and  Taylor,  1995).  Numerical  experiments  by  McLaughlin  and  Jih  (1986, 
1987)  and  Jih  (1993b,  1995)  readily  demonstrate  that  rough  topography  and  random  heterogeneity  can 
scatter  significant  Rg  energy  into  body  waves.  An  additional  Rg  -io-Lg  conversion  mechanism  has  just 
been  reported  (Jih,  1 995).  This  new  Lg  excitation  mechanism  is  rather  interesting  and  worth  some  dis¬ 
cussion  here.  It  might  be  anticipated  intuitively  that  the  net  effect  of  anelasticity  in  the  surface  layers  is 
solely  to  reduce  the  amplitude  of  incident  seismic  waves.  The  Rg  wave  is  particularly  susceptible  to 
such  a  mechanism  since  it  is  confined  in  the  uppermost  crustal  layers  where  the  anelastic  attenuation  is 
often  very  strong.  However,  if  the  anelastic  attenuating  layer  is  only  ftick  enough  to  dissipate  the  retro¬ 
grade  rolling  near  the  surface,  then  the  free  surface  would  behave  asymptotically  like  a  fixed  point. 
Beyond  a  certain  distance,  the  fundamental  mode  can  no  longer  be  sustained  by  such  a  waveguide,  and 
accordingly  any  undissipated  Rg  energy  would  have  to  propagate  in  other  wave  types  or  modes.  This  is 
very  similar  to  the  situation  of  a  solid  half  space  with  a  rigid  surface.  There  is  no  corresponding  funda¬ 
mental  mode  Rayleigh  wave  in  such  a  structure  (Aid  and  Richards,  1980,  pp.189).  Jih  (1995)  shows 
that  this  process  couples  the  undissipated  Rg  energy  into  pure  shear  waves  or  higher  modes,  depend¬ 
ing  on  the  complexity  of  the  structure.  In  terms  of  the  Rg  -to-S  to  Rg  -to-P  ratio,  this  process  appears 
to  be  far  more  efficient  than  other  near-surface  Rg  scattering  mechanisms.  Furthermore,  in  this  pro¬ 
ems,  the  Rg  spectrum  is  naturally  imjxinted  onto  the  converted  S  waves,  which  could  help  to  explain 
some  recently  observed  spectral  characteristics  of  Lg  waves. 

Two  LFD  experiments  have  been  conducted  here  to  illustrate  these  Rg  -related  findings.  The  first 
model  is  an  anelastic  half  space  with  compressional  and  shear  velocities  of  5.02  and  2.898  km/sec, 
respectively.  The  damping  factor  is  set  to  0.9995  at  every  grid  point.  The  vertical-component  snapshots 
of  Rg  wave  propagation  in  this  anelastic  half  space  are  shown  in  Figure  11.  Figure  12  shows  vertical- 
component  record  sections  of  planar  Rg  waves  recorded  at  surface  sensors  (top)  and  5-km  depth  (bot¬ 
tom)  for  this  anelastic  half  space  model.  The  Rg  energy  is  confined  in  the  near-surface  layer  and  very 
little  energy  is  detected  at  a  depth  of  5  km,  as  expected.  The  fastest  Rg  phase  has  a  group  velocity  of 
2.33  km/sec,  which  is  slower  than  that  of  Rg  in  the  elastic  case  (Figure  12).  The  slowest  Rg  phase  has 
a  velocity  of  only  about  1  km/sec  (Figure  11).  Figures  13  and  14  show  the  amplitude  spectra  of  traces 
No.  5  and  No.  20,  respectively.  These  two  surface  sensors  are  15  km  apart.  The  spectral  ratio  is  very 
flat  over  the  frequency  band  between  2  and  10  Hz  (Figure  15). 

In  the  second  test,  all  the  parameters  are  identical  to  the  first  one  except  that  the  low-Q  layer  is 
limited  to  the  uppermost  500  meters  only  (see  Rgure  16).  Due  to  this  peculiar  structure,  some  Rg 
energy  starts  to  detach  from  the  Rayleigh  mode  and  eventually  is  converted  into  primarily  pure  shear 
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waves.  These  shear  waves  propagate  at  a  F>ost-crrtical  angle  and  hence  they  are  excellent  candidates 
for  regional  shear  phases  such  as  Sg ,  Lg  or  Lg  coda.  Figure  1 7  shows  the  vertical-component  record 
sections  recorded  at  the  free  surface  (top)  and  at  5-km  depth  (bottom)  for  this  hybrid  model.  The  con¬ 
verted  shear  waves  (recorded  at  a  depth  of  5  km)  have  peak  amplitudes  actually  greater  than  those  of 
Rayleigh  waves  recorded  at  the  surface  sensors.  It  is  also  interesting  to  note  that,  while  the  Rg  wave 
train  travels  at  a  velocity  of  2.58  km/sec  along  the  free  surface,  the  converted  shear  waves  have  a  velo¬ 
city  about  3.3  km/sec.  This  provides  yet  another  piece  of  evidence  that  the  converted  shear  waves  are 
good  candidates  for  the  Lg  -type  of  phase.  Figures  18  through  20  show  the  spectra  and  ratio  of  traces 
No.  5  and  No.  20  for  this  hybrid  model.  A  synthetic  seismogram  recorded  at  5-km  depth  is  compared 
against  the  reference  surface  sensor  No.  5  (Rgures  21  and  22).  All  the  spectral  analysis  and  Q  results 
are  summarized  in  Table  5  below.  Even  though  the  same  damping  factor  of  0.9995  was  used  in  these 
two  experiments,  the  resulting  Q  values,  and  the  shape  of  transmitted  Rg  waves,  as  well  as  the  wave 
field,  all  look  very  different.  Apparently  the  thickness  of  the  attenuating  layer  is  very  important.  This  is 
particularly  true  for  Rg  waves. 

However,  it  should  be  emphasized  that  the  role  of  anelastic  attenuation  in  this  new  Rg  -to-SVL^ 
conversion  mechanism  is  to  alter  the  eigenfunction  corresponding  to  this  peculiar  structure,  through  dis¬ 
sipation  of  Rg  energy  near  the  free  surface.  The  low  Q  layer  itself  does  not  preserve  energy,  nor  does 
it  re-direct  any  energy  to  escape  from  the  layer.  It  is  the  energy  in  the  Rg  root,  which  was  origindly  rol¬ 
ling  underneath  the  top  layer  and  was  never  dissipated  by  the  absorption  mechanism  on  the  top,  that 
eventually  becomes  detached  from  the  energy  trapped  in  the  top  layer.  The  energy  in  the  Rg  root  is 
preserved  and  converted  to  S  or  .  Therefore,  this  Rg  -to-S  or  Rg  -to-Lg  mechanism  should  be 
regarded  as  a  mixture  of  absorption  and  scattering. 
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Table  5.  Rg  Propagation  in  Anelastic  Crustal  Models 


Damper 

Sensor  Location  (km) 

Results  at  2  Hz 

Results  at  4 

Hz 

G,  Thickness  (km) 

Range 

Depth 

A/Ao 

Q(2Hz) 

y  (2  Hz) 

A/Ao 

Q(4  Hz) 

Y(4Hz) 

0.9995, 

9 

0 

0.729 

63 

0.0420 

0.729 

150 

0.0351 

0.9995,  '5° 

10 

0 

0.691 

62 

0.0422 

0.692 

143 

0.0368 

0.9995, 

11 

0 

0.655 

60 

0.0438 

0.656 

138 

0.0383 

0.9995, 

12 

0 

0.620 

56 

0.0469 

0.621 

133 

0.0397 

0.9995,  oo 

13 

0 

0.587 

56 

0.0469 

0.588 

129 

0.0408 

0.9995, 

14 

0 

0.557 

57 

0.0465 

0.559 

127 

0.0416 

0.9995, 

15 

0 

0.527 

54 

0.0485 

0.529 

124 

0.0424 

0.9995, 

16 

0 

0.498 

53 

0.0500 

0.501 

122 

0.0432 

0.9995, 

17 

0 

0.473 

54 

0.0486 

0.476 

121 

0.0437 

0.9995, 

18 

0 

0.454 

54 

0.0487 

0.457 

121 

0.0435 

0.9995,  0.5 

9 

0 

0.543 

26 

0.0931 

0.555 

73 

0.0653 

0.9995,  0.5 

10 

0 

0.513 

26 

0.0920 

0.526 

74 

0.0642 

0.9995,  0.5 

11 

0 

0.487 

26 

0.0903 

0.501 

76 

0.0628 

0.9995,  0.5 

12 

0 

0.464 

27 

0.0885 

0.479 

78 

0.0614 

0.9995,  0.5 

13 

0 

0.444 

28 

0.0865 

0.459 

80 

0.0598 

0.9995,  0.5 

14 

0 

0.426 

28 

0.0846 

0.442 

82 

0.0584 

0.9995,  0.5 

15 

0 

0.410 

29 

0.0826 

0.425 

84 

0.0570 

0.9995,  0.5 

16 

0 

0.395 

30 

0.0806 

0.411 

86 

0.0556 

0.9995,  0.5 

17 

0 

0.381 

30 

0.0785 

0.397 

88 

0.0544 

0.9995,  0.5 

18 

0 

0.366 

31 

0.0770 

0.382 

89 

0.9995,  0.5 

9 

5 

0.243 

14 

0.1341 

0.239 

23 

0.1590 

0.9995,  0.5 

10 

5 

0.282 

17 

0.1060 

0.277 

29 

0.1283 

0.9995,  0.5 

11 

5 

0.311 

21 

0.0879 

0.306 

34 

0.1077 

0.9995,  0.5 

12 

5 

0.330 

24 

0.0761 

0.324 

39 

0.0938 

0.9995,  0.5 

13 

5 

0.340 

27 

0.0682 

0.335 

44 

0.0842 

0.9995,  0.5 

14 

5 

0.346 

29 

0.0626 

0.341 

48 

0.0770 

0.9995,  0.5 

15 

5 

0.348 

32 

0.0585 

0.342 

52 

0.0715 

0.9995,  0.5 

16 

5 

0.346 

33 

0.0554 

0.341 

55 

0.0673 

0.9995,  0.5 

17 

5 

0.343 

35 

0.0529 

0.338 

58 

0.0638 

0.9995,  0.5 

18 

5 

0.338 

36 

0.0511 

0.333 

60 

0.0611 
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FDstep  0 

VERT  0.000S9C,  1  800  1  200 

XTIcMarfc=  1  km;  Z  Tic  Mark  =  1km 
-OJBSE’M,  0.19Ef01 


FDstop  400 

VERT  1.000sec,  1  800  1  200 

XTIcMarks  1  km;  Z  Tic  Mark  =  1km 
-0.11E+01,  0.14E+01 


FDatap  800 

VERT  2J00099C,  1  800  1  200 

XTIcMark=  1  km;  Z Tic Maik  =  1km 
-0.11  E+01,  0B7E*O0 


FDstep  1200 

VERT  SXKXteac,  1  800  1  200 
XTIcMark=  1  km;  Z  Tic  Mark  =  1km 
-O^SE-i^,  0.80E+00 


FDstep  1800 

VERT  4.S00sec,  1  800  1  200 

XTlcMark=  1  km;  Z  Tic  Mark  =  1km 
-0.72E<^00,  0.62E<^00 


FDstop  2400 

VERT  6.000sac,  1  800  1  200 

XTIcMaiks  1  km;  ZTIcMarks  1km 
-0.55E-^00,  0.51  Et^ 


FDstep  3000 

VERT  7.500SOC,  1  800  1  200 

XTlcMark=  1  km;  Z  Tic  Mark  =  1km 
-0>43Ei^,  0.41E-f00 


LFD  Simulation  of  Rg  Propagation:  Model  Q.9995,600 

Figure  11.  The  vertical-component  snapshots  of  Rg  wave  propagation  in  an  anelastic  half  space.  The 
compressional  and  shear  velocities  are  5.02  and  2.898  km/sec,  respectively.  The  damping  factor  is  0.9995 
at  every  grid  point  The  intrinsic  attenuation  in  the  half  space  causes  a  severe  dispersion,  which  has  been 
shown  to  be  a  necessary  condition  of  causality. 
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Rg+HS+Q;  G(HS)=0.9995,  Z  component,  sensor  depth  =  5km 


Figure  12.  Vertical-component  record  sections  of  planar  Rg  waves  recorded  at  surface  sensors  (top)  and  5-km  depth 
(bottom)  for  the  anelastic  half  space  model.  The  Rg  energy  is  confined  in  the  near-surface  layer  and  very  little  energy  is 
detected  at  a  depth  of  5  km. 
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Specalc+Waveshaper 
version  9502,  JB+TWM+RJ 

Window  1 024 
Signal  1  1024 

10.00%  Cosine  Taper 
Taper  Both  Ends 

41  R  Smoothing 
No  Decimation 

Deconvolution  Flag=0.  - 

Waveshaper  Fiag=0. 

Reject:  <  O.OOHz 
Reject:  >1 1  .OOHz 
Pass  Band=  0.01,10.00Hz 

Final  Taper=1 0.00%, 10.00% 

Model  Q.9995,600 
VERT.5.  Rg 


1 .28  sec  1 .341 0  Input  Trace  1 .28  sec  0.0013  Output  T race 


Rgure  13.  Amplitude  spectrum  of  the  vertical-component  displacement  recorded  at  sensor  No.  5  of  the  anelastic  half¬ 
space  model  with  a  damping  factor  of  0.9995.  The  sensor  is  located  near  the  left  portion  of  the  grid. 
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Specalc+Waveshaper 

version  9502,  JB+TWM+RJ 

Window  1024 
Signal  2641  1024 

1 0.00%  Cosine  Taper 
Taper  Both  Ends 

41  Pt  Smoothing 
No  Decimation 

Deconvolution  Flag=0. 

Waveshaper  Flag=0. 

Reject:  <  O.OOHz 
Reject:  >1 1  .OOHz 
Pass  Band=  0.01 ,10.00Hz 

Final  Taper=1 0.00%, 1 0.00% 

Model  Q.9995,600 
VERT.20,  Rg 


1 .28  sec  0.4050  Input  T race 


I _ l 

1 .28  sec  0.0004  Output  T  race 


Figure  14.  Amplitude  spectrum  of  the  vertical-component  displacement  recorded  at  surface  sensor  No.  20  of  the  ane- 
lastic  half-space  model.  The  sensor  is  located  near  the  right  portion  of  the  grid. 
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TGAL  RATIO  vers  1 .5  (JB+TWM+RAW+RJ) 

Fri  Nov  10  09:44:02  1995 

Model  Q.9995,600,  ratio:  No.20/  No.5 

NOISE  POWER  NOT  SUBTRACTED  FROM  SIGNAL  POWER 

S/N  POWER  THRESHOLD  =  2.0 

7  POINT  SMOOTHING 


O 

o 


FREQ  BAND  SLOPE  +/-  Y INTERC  AV INTEG  AV  RATIO  SD 
2.00-4.00  6.314e-03  2.437e-03  4.340e-01  -3.435e-01  -3.427e-01  5.153e-03 

4.00-6.00  1.000e-06  1.245e-06  4.570e-01  -3.401  e-01  -3.401  e-01  O.OOOe+00 

6.00-8.00  -8.350e-05  1.612e-05  4.575e-01  -3.402e-01  -3.402e-01  O.OOOe+00 

3.00-8.00  -4.565e-06  8.843e-06  4.570e-01  -3.402e-01  -3.402e-01  1.410e-04 

2.00-10.00  5.078e-04  2.300e-04  4.531  e-01  -3.408e-01  -3.407e-01  2.651  e-03 


Figure  15.  Rg  amplitude  ratio  of  synthetic  seismograms  No.  20  to  No.  5  of  the  anelastic  half-space  mode! 
with  a  damping  factor  0.9995.  These  two  sensors  are  15  km  apart  (see  Figure  12). 
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FDstep  0 

VERT  0.00088C,  1  800  1  200 

XTlcMarks  1  km;  Z Tic  Marie  =  1km 
•OJSaEM,  0.19E401 


FDatap  400 

VERT  IXlOOsac,  1  800  1  200 
XTIcMaric=:  1  km;  ZTIcMarkr  1km 
-0.11  E+01,  0.66E^ 


FDstap  800 

VERT  2J)00sec,  1  800  1  200 

XTIcMaricr  1  km;  Z Tic  Mark  =  1km 
-0.65E400,  0.74E400 


FDstop  1200 

VERT  3UX)Omc,  1  800  1  200 
XTIcMaricr  1  km;  Z  Tic  Marie  r  1km 
-0.55E-f00,  0.57E+00 


FDstep  1800 

VERT  43003ec,  1  800  1  200 

XTIcMaricr  1  km;  ZTIcMarkr  1  km 
’OA2E*0O,  OMEM 


FDstep  2400 

VERT  6.0008ec,  1  800  1  200 

XTIcMaricr  1  km;  Z  Tic  Mark  r  1km 
-0.36E<i-00,  0.32E+00 


FDstep  3000 

VERT  7.50asec,  1  800  1  200 

XTIcMaricr  1  km;  ZTIcMarkr  1km 
-0.11  EtDO,  0.22E+00 


LFD  Simulation  of  Rg  Propagation:  Model  Q.9995,22 

Figure  16.  The  vertical-component  snapshots  of  Rg  wave  propagation  in  an  elastic  half  space  with  an 
anelastic  layer  of  500  meters  thick  lying  over  it.  Due  to  this  peculiar  structure,  some  Rg  energy  appears  to 
have  gradually  detached  from  Rayleigh  mode  and  been  converted  into  primarily  pure  shear  waves. 
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Rg+HS+Q;  G(top  500m)=0.9995.  Z  component,  sensor  depth  =  5km 


Figure  17.  Vertical-component  record  sections  of  planar  Rg  waves  recorded  at  surface  sensors  (top)  and  5-km  depth 
(bottom)  for  the  elastic  half  space  with  a  shallow  anelastic  layer  lying  over  it  The  converted  shear  waves  have  peak 
amplitudes  greater  than  those  of  Rayleigh  waves  recorded  at  the  surface  sensors. 
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Specalc+Waveshaper 

version  9502,  JB+TWM+RJ 

Window  1 024 
Signal  1  1024 

10.00%  Cosine  Taper 
Taper  Both  Ends 

41  R  Smoothing 
No  Decimation 

Deconvolution  Fiag=0. 

Waveshaper  Flag=0. 

Reject:  <  O.OOHz 
Reject:  >1 1  .OOHz 
Pass  Band=  0.01,10.00Hz 

Final  Taper=1 0.00%, 10.00% 

Model  Q.9995,22 
VERT.5,  Rg 


1 .28  sec  1 .0950  Input  Trace  1 .28  sec  0.001 1  Output  T race 


Rgure  18.  Amplitude  spectrum  of  the  vertical-component  displacement  recorded  at  sensor  No.  5  of  the  anelastic  model. 
This  surface  sensor  is  located  near  the  left  portion  of  the  grid. 
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Specalc+Waveshaper 
version  9502,  JB+TWM+RJ 

Window  1 024 
Signal  2191  1024 

1 0.00%  Cosine  Taper 
Taper  Both  Ends 

41  R  Smoothing 
No  Decimation 

Deconvolution  Flag^O. 

Waveshaper  Flag=0. 

Reject:  <  O.OOHz 
Reject:  >1 1  .OOHz 
Pass  Band=  0.01 ,10.00Hz 

Final  Taper=1 0.00%, 10.00% 

Model  Q.9995,22 
VERT.20,  Rg 


I _ I 

1 .28  sec  0.2938  Input  Trace 


I _ I 

1 .28  sec  0.0003  Output  T race 


Rgure  19.  Amplitude  spectrum  of  the  vertical-component  displacement  recorded  at  sensor  No.  20  of  the  half-space 
model  with  a  shallow  attenuating  layer  lying  over  it.  This  surface  sensor  is  located  near  the  right  portion  of  the  grid. 
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TGAL  RATIO  vers  1.5  (JB+TWM+RAW+RJ) 

Fri  Nov  10  09:42:16  1995 

Model  Q.9995,22.  ratio:  No.20  /  No.5 

NOISE  POWER  NOT  SUBTRACTED  FROM  SIGNAL  POWER 
S/N  POWER  THRESHOLD  =  2.0 
7  POINT  SMOOTHING 


O 


FREQUENCY  (Hz) 

FREQ  BAND  SLOPE  +/-  Y INTERC  AV INTEG  AV  RATIO  SD 
2.00-4.00  4.366e-02  1.601e-02  2.676e-01  -4.415e-01  -4.361  e-01  3.487e-02 
4.00-6.00  1.500e-06  1.327e-06  3.821  e-01  -4.179e-01  -4.179e-01  NaN 
6.00-8.00  1.800e-05  7.278e-06  3.820e-01  -4.178e-01  -4.178e-01  1.726e-04 
3.00-8.00  5.249e-05  2.468e-05  3.818e-01  -4.179e-01  -4.179e-01  1.410e-04 
2.00-10.00  3.675e-03  1.551e-03  3.597e-01  -4.220e-01  -4.218e-01  1.814e-02 


Rgure  20.  Rg  amplitude  ratio  of  synthetic  seismograms  No.  20  to  No.  5  of  the  anelastic  half-space 
model.  These  two  sensors  are  15  km  apart  (see  Figure  17). 
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Specalc+Waveshaper 

version  9502.  JB+TWM+RJ 

Window  1 024 
Signal  1756  1024 

1 0.00%  Cosine  Taper 
Taper  Both  Ends 

41  R  Smoothing 
No  Decimation 

Deconvolution  Flag=0. 

Waveshaper  Flag=0. 

Reject:  <  O.OOHz 
Reject:  >1 1  .OOHz 
Pass  Band=  0.01,10.00Hz 

Final  Taper=1 0.00%, 10.00% 

Model  Q.9995,22 
VERT.40,  Rg 


I _ I 

1.28  sec  0.3371  Input  Trace 


1 .28  sec  0.0003  Output  T race 


Figure  21 .  Amplitude  spectrum  of  the  vertical-component  displacement  recorded  at  the  sensor  No.  40  of  the  half-space 
model  with  a  shallow  attenuating  layer  lying  over  it.  This  sensor  is  located  near  the  right  portion  of  the  grid,  at  a  depth  of 
5  km. 
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TGAL  RATIO  vers  1 .5  (JB+TWM+RAW+RJ) 

Fri  Nov  10  09:43:09  1995 

Model  Q.9995,22,  ratio:  No.40  /  No.5 

NOISE  POWER  NOT  SUBTRACTED  FROM  SIGNAL  POWER 
S/N  POWER  THRESHOLD  =  2.0 
7  POINT  SMOOTHING 


O 


FREQ  BAND  SLOPE  +/-  Y INTERC  AV INTEG  AV  RATIO  SD 
2.00-4.00  -1.563e-02  5.172e-03  3.785e-01  -4.688e-01  -4.708e-01  1.201e-02 

4.00-6.00  1.500e-06  2.869e-06  3.330e-01  -4.775e-01  -4.775e-01  O.OOOe+00 

6.00-8.00  -6.570e-04  2.051  e-04  3.364e-01  -4.778e-01  -4.778e-01  4.883e-04 

3.00-8.00  -1.600e-04  4.867e-05  3.336e-01  -4.776e-01  -4.776e-01  3.729e-04 

2.00-10.00  -3.707e-03  6.229e-04  3.488e-01  -4.797e-01  -4.799e-01  1.072e-02 


Rgure  22.  Rg  amplitude  ratio  of  synthetic  seismograms  No.  40  to  No.  5  of  the  half-space  model  with  a 
shallow  attenuating  layer  lying  over  it.  Trace  No.  40  is  15  km  away  from  the  surface  sensor  No.  5.  At  a 
depth  of  5  tan.  sensor  No.  40  only  records  the  converted  shear  waves  that  travel  at  a  postcritical  angle. 
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6.  DISCUSSION  AND  CONCLUSION 


A  simple,  causal  method  for  incorporating  anelastic  attenuation  into  the  LFD  method  is  presented 
in  this  report  The  method  is  simple  to  implement.  Testing  of  this  algorithm  with  various  seismic  waves 
demonstrates  that  this  damping  operator  provides  the  expected  intrinsic  attenuation  over  a  rather  broad 
frequency  band.  This  ad  hoc  damping  gives  an  intrinsic  attenuation  with  the  quality  factor  inaeasing 
linearly  with  the  frequency.  Also,  in  the  P-SV  situation,  the  resulting  Op  and  Os  are  approximately  the 
same.  If  a  frequency  dependence  other  than  the  linear  one  is  desired,  then  several  separate  finite- 
difference  simulations  need  to  be  carried  out  for  each  individual  frequency-0  pair  of  interest.  The 
shortcomings  of  this  procedure  are  outweighed  by  the  simplicity.  More  importantly,  this  procedure 
preserves  the  causality.  Using  the  algorithm  described  in  this  study,  it  is  easier  to  relate  and  equate  the 
propagation  effects  due  to  small-scale  random  heterogeneities  (and  other  large-scale  structural  varia¬ 
tions)  with  those  due  to  the  anelastic  attenuation  (see  Jih,  1995, 1996). 

The  procedure  described  in  this  study  is  based  on  that  of  Cerjan  et  at.  (1985).  It  is  shown  that 
this  simple  procedure  results  in  an  amplitude  that  decays  exponentially  with  traveling  distance,  which 
happens  to  be  an  approximate  solution  to  the  telegraphy  equation  (Levander,  1985a).  The  damping 
effect  added  to  the  wave  equation  is  equivalent  to,  at  least  asymptotically,  the  friction  term  Levander 
(1985a)  added  to  the  wave  equation.  Thus  the  empirical  damping  also  has  a  theoretical  justification. 
Since  our  damping  term,  G,  is  meant  to  represent  the  spatial  attenuation,  which  is  a  very  localized 
phenomenon,  the  damping  term  could  vary  from  grid  to  grid.  However,  if  it  is  used  as  a  supplementary 
absorbing  boundary  condition,  there  may  be  a  gradient  in  G  that  would  provide  a  better  effect  (see 
Sochaki  et  at.,  1987). 

All  the  examples  tested  in  this  study  consistently  show  that,  the  performance  of  this  damping  pro¬ 
cedure  as  an  attenuation  operator  degrades  somewhat  at  very  low  frequencies.  There  could  be  damp¬ 
ing  or  amplification  of  the  signal  at  the  low  frequency  end,  a  phenomenon  which  Levander  (1985^  also 
observed  in  utilizing  the  telegraphy  equation  as  an  auxiliary  absorbing  boundary.  Nevertheless,  as 
Levander  (1985a)  pointed  out  the  introduction  of  a  low-frequency  reflection  (or  amplification)  is  more 
desirable  than  a  high-frequency  reflection  (or  amplification).  It  is  appropriate  to  regard  this  damping  pro¬ 
cedure  (as  well  as  Levanderis  use  of  the  telegraphy  equation)  as  an  ad  hoc  measure,  useful  until  a 
more  robust  attenuation  operator  or  absorbing  boundary  condition  is  developed. 
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